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ABSTRACT 


An   investigation   was   conducted  to  determine  the 

feasibility  of  obtaining  the  short  period  and   phugcid 

stability     derivatives     from    one    maneuver, 

"simultaneously".  It  was  concluded   that   the   maximum 

likelihood   identification   program   SCIDNT   I   showed 

great   promise   in   obtaining   the   short   period   and 

phugoid     stability     derivatives     from     one 

maneuver.  Extraction  of   the   short   period  stability 

parameters   in    the   presence   of   the   phugoid   was 

easy,  straight  forward,  and  yielded  results  similar  to 

those  obtained  from  pure  short  period  data.  Estimation 

of  the  phugoid  stability  parameters  was  possible   when 

they     were    estimated     in    conjunction     with 

Z  ,  X  ,  H  ,  and  9  .  It  was  recommended  that  a  new   set 

o    o    o       o 
of   data   be   obtained  at  a  flight  condition  where  the 

phugoid  is  at  least  moderately  damped   and   that   this 

data  be  analyzed  to  resolve  the  present  anomalies. 
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I.   INTRODUCTION 


A.   BACKGROUND 

Flight  testing  of  military  aircraft  has  always  been  a 
time  consuming  and  expensive  proposition.  The  relatively 
recent  introduction  of  the  new  Flying  Qualities 
Specification  (MIL-F-S785B)  has  introduced  many  new 
parametric  requirements  which  require  that  many  of  the 
aircraft  airframe  parameters  be  determined  to  insure 
specification  compliance.  As  a  result  of  these  new 
requirements  the  Naval  Air  Test  Center  (NATC)  initiated  a 
research  program  to  develop  airframe  parameter 
identification  technology  for  use  in  flight  testing  Navy 
aircraft.  This  program  was  initiated  in  1971  between  the 
Naval  Air  Systems  Command  (NAVAIR)  and  NATC.  At  the  same 
time  the  Office  of  Naval  Research  (ONR)  issued  a  contract  to 
Systems  Control,  Inc.  (SCI)  to  advance  the  state-of-the-art 
in  parameter  identification.  In  1973  a  concentrated  effort 
was  begun  by  the  Navy  to  continue  the  development  of  this 
new  technology  by  forming  a  joint  program  between 
NATC,  NAVAIR,  ONR  and  SCI.  As  a  result  of  this  development  a 
digital  FORTRAN  computer  program,  SCIDNT-I  (hereafter 
refered  to  as  SCIDNT)  was  written  for  NATC  to  extract 
aircraft  stability  and  control  coefficients  from  flight  test 
data,  Ref.  1.  it  was  hoped  that  by  using  this  computer 
program  it  would  be  possible  to  significantly  reduce  the 
flight  test  time  in  the  air,  and  the  ground  data  reduction. 
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B.   PURPOSE 

The  purpose  of  this  research  was  to  develop  parameter 
identification  test  methods  for  determining  aircraft  phugoid 
and  short  period  characteristics  from  one  aircraft 
maneuver.  This  research  further  entailed  the  determination 
of  proper  aircraft  maneuvers/pilot  inputs  and  detailed  data 
analysis  proceedures  using  the  advanced  maximum  likelihood 
identification  technigue. 


C.   DESCRIPTION  OF  AIRPLANE 

The  test  airplane  was  demonstration  F-14A  ship  number  8 
BU.  NO.  157S87.  The  airplane  was  essentially  typical  of  a 
production  F-14A  airplane  relative  to  these  tests.  An 
aircraft  description  is  contained  in  Appendix  A.  Photographs 
and  a  three  view  drawing  of  the  test  airplane  are  shown  in 
Appendix  A  ,  Figures  1,  2  and  3. 


D.   tiAXIMUM  LIKELIHOOD  PARAMETER  IDENTIFICATION  THEORY 

A  widely  used  rule  for  choosing  the  sample  statistic  to 
be  used  in  parameter  estimation  is  called  the  method  of 
maximum  likelihood  Ref.  2.  This  method,  in  addition  to 
prescribing  the  statistic  which  should  be  used,  provides  an 
approximation  of  the  distribution  of  that  statistic  so  that 
approximate  confidence  intervals  may  be  constructed.  The 
method  of  maximum  likelihood  makes  use  of  the 
sample- likelihood  function.  Assume  we  observe  or  record  a 
set  of  observations  on  a  system,  whose  model  has  p  unknown 
parameters  6.  -We  then  can  write  the  joint  probability 
distribution  of  a  random  sample  for  any  given  set  of   values 
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of  the  parameters  6  from  a  feasible  set  8.  We  then  can 
assign  a  probability  p(z|6)  to  each  outcome  z.  If  we 
observe  or  measure  the  outcome  of  an  actual  experiment  to  be 
z  we  then  would  like  to  know  which  set  or  sets  of  values  of 
6  might  have  resulted  in  these  observations.  This  concept  is 
contained  in  the  likelihood  function  ac  (9|z).  This  function 
is  of  fundamental  importance  in  estimation  theory  because  of 
the  likelihood  principle  of  Fisher  and  others  which  states 
that  if  the  system  model  is  correct,  all  information  about 
unknown  parameters  is  contained  in  the  likelihood 
function.  The  maximum  likelihood  estimator  of  9  is  the  value 
6  which  causes  the  likelihood  function  aC  (9|z)  to  be  a 
maximum  [2J. 


%   =  max   £  (6|z) 
B€9 


From  a   purely   mechanical   viewpoint   this   method   can   be 
summarized  as  follows: 


Find  the  probability  density  functions  of  the 
observations  or  measured  values  for  all  possible 
combinations  of  the  unjcnown  parameter 
values.  Select  the  density  function  whose  value  is 
highest  among  all  density  functions  at  the 
measured  values  of  the  observations.  The 
corresponding  parameter  values  are  the  maximum 
likelihood  estimates. 


Assume  9  can  take   three   possible   values:   9,9,  and 

12 

9  .  Further  assume  that  the  probability  density  functions  of 

observations  z  for  these  three  values  of  9  are  as   shown   in 
Figure   1.  Then,  if   the   actual  observation  is  z,  9   is  the 

2 

maximum  likelihood  estimate  of  9. 
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p(z|6) 


observed 
Figure  1.   Maximum  Likelihood  Estimates 


Owing  to  the  monotonic,  one-to-one  relationship  between 
the  likelihood  function  and  its  logarithm,  i.e.  they  have  a 
maximum  at  the  same  value  8,  it  is  more  convient  tc  work 
with  the  negative  of  the  logarithm  of  the  likelihood 
function.  Starting  from  a  priori  values,  the  parameters  are 
updated  so  that  the  value  of  the  negative  log- likelihood 
function  J  at  the  observed  value  of  the  outputs  decreases 
continuously . 

As  sample  statistics  or  functions  of  the 
observations,  the  estimators  are,  before  the  experiment, 
random  variables  and  can  be  studied  as  such.  In  fact,  the 
major  advantage  of  maximum  likelihood  estimators  is  that 
their  properties  have  been  thoroughly  studied  and  are  well 
known,  virtually  independently  of  the  particular  density 
function  under  consideration.  Additionally,  the  maximum 
likelihood  method  can  be  used  with  linear  or  nonlinear 
models  in  the  presence  of  process  and  measurement  noise.  The 
estimator  of  6  is   approximately   normally   distributed   for 
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large     n,  and    is   asymptotically    unbiased,  efficient, 
consistent,  and  sufficient,  [2], 

The   details   of   the   maximum  likelihood  identification 
proceedure  are  presented  in  the  following  section. 


E.   MAXIMUM  LIKELIHOOD  IDENTIFICATION  OP  PARAMETERS 

The  derivation  follows  closely  that  presented  in  Ref. 
3,  and  is  presented  herein  for  completeness  and  reader 
convenience.  The  initial  equations  will  be  concerned  with 
general  ncnlinear  systems. 

Consider  the  general  nonlinear  aircraft  equations  of 
motion 


x  =  f(x,  u,  6,  t)  +  P(6,  t)w  0  <  t  <  T        (  1  ) 


E(x(0))  =  x  (6) 
o 


E{(x(0)  -  x  (6))  (x{0)  -  x  (0))  }  =  P  (9)  (  2  ) 

o  o  o 

where 


x(t)  is  n  x  1  state  vector 

u(t)  is  1  x  1  input  vector 

6  is  p  x  1  vector  of  unknown  parameters 

P(9,t)  is  n  x  q  process  noise  distribution  matrix 

w  (t)  is  q  x  1  random  process  noise  vector 


Sets  of  m  measurements  y (t  )  are  taken  at  discrete  times  t 

k  k 
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y(t  )  =  h(x(t  ),  u(t  ),  e,  t  )  +  v(t  )        (  3  ) 

k        k      k       k       k 

k  =  1,  2,  3, tN 


w(t)   and   v  (t  )   are  Gaussian   random  noises  [2],  with  the 
k 


following  properties 


E(w(t))  =  0      E(v(t  ))  =  0 

k 

T 
E(w(t)  v  (t  ))  =0 
k 


T 
E(w(t)  w  (T))  =  Q(6,t)  6   (t-T) 


2(v(t  .)  v  (t  ))  =  R(9,t  )    O  (  4  ) 

3       *  j    jk 

The   unknown   parameters   are  assumed   to  occur  in  the 

functions  f  and  h  and  in   matrices  P,  Q,  R,  P  ,    and   x  .  In 

1  o        o 

the  following  analysis,  the  model  and  the  functional  form  of 

f  and  h  is  assumed  known  correctly. 

The     set    of   observations    y  (t  ),y(t  ),....,  y  (t  ) 

1        2  N 

constitutes  the   outcome   Z   in   this   case.  The   likelihood 

function  for  8,  which  has  the  same  form  as  the  probability 
of  the  outcome  z  for  a  certain  value  of  parameters  6,  is 
given  by 


g   (6|z)  =  p(z|9) 


-  P<Y<t  )  ,  y(t  )  ,...,y(t  )  |e) 

12  n 
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^  (G|z)  =  p(Y  |6)  (  5  ) 

N 


where 


i  =  £y(t  ) ,...,y  (t  )j ,  k  =  1,  2,  ...,n 

k       l         k 


p(Y  |0)  =  p(y(t  )  |Y    ,6).  p(Y    |9) 
N  N    N-1  N-1 


=  p^'V'Vre)  p(y(tN-i,,YN-2'  e)  P<V2,9) 


which  by  successive  application  of  Bayes  rule,  [2],  becomes 


p(Y  |9)  =  TT  p(y(t.)  |Y.   ,  e)  (  6  ) 

N        '  '        1    1-1 
i=1 


The  log- likelihood  function  is 


N 

log  (&(e\z))    =    )        log  {p(y(t  )  |Y    ,  6)  }  +  constant   (  7  ) 

*-•  i        i-1 

i-1 


To    find   the    probability   distribution   of    y(t    )    given    Y  and 

i         i-1 

6,  the  mean  value  and  covariance  are  determined  first. 


E(y(t.)  iy.   ,  e)  A  y(i|i-i)  (  8  ) 

i   1-1 


The  expected  value  or  the  mean  is  the  best  possible  estimate 
of   measurements  at  a  point  given  the  measurements  up  to  and 
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including  the  previous  point. 


cov(y(t.)|Y.   ,e)  =  E{(y(t.)  -  y  (i|i-1J )  (y  (t .)  -  y(i|i-1))  } 

1   i-1  1  i 


A  E{V(i)  \T(i)    } 


6  B(i)  (  9  ) 


V*'  (i)  are  called  the  innovations  at  point  i  and  B(i)  is   the 
innovations  covariance.  Since 


y  (t.)  -  E(y(t.)  |Y.   ,9)  =  V  (i)  (  10  ) 

l  ii-1 


it   follows   that   y (t  )   given   Y      and   6   have  the  same 

i  i-1 

distribution   as  v  (i).  It   has   been   shown   that   as   the 

sampling    rate   is   increased,  the   innovations  v   (i)   tend 

towards  having  a  Gaussian  density.  Assuming   a   sufficiently 

high   sampling    rate,  the    distribution   of  V (i)    and, 

therefore,  y(t  )  given  Y     and  6  is  Gaussian,  i.e., 
i         i-1 


exp 


p(y(t.)  iy.  ,9)  =• 


|-— -  V(i)T  B-i(i)  ^(i)  ! 


i    i-1  m/2       1/2 

(27T)     |B(i)| 


(  11  ) 
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log  {p(y(t.)|Y.   ,6)}  =  ~  —  V  (i)  B-»(i)  ^(i) 
11-1  2 


1 

—  log|B(i)|  +  constant 
2 


(  12  ) 


The  log-likelihood  function  of  equation  (  7  )  can  be  written 
as 


1  r-     T 

log  (  s£(6|z))  =  -  —  2^  K  (i)B-Mi)V(i) 


2 
i=1 


♦  log|B(i)  |}  (  13  ) 


An  estimate  of  the  unknown  parameters  is  obtained  by 
maximizing  the  likelihood  function  or  the  log-likelihood 
function  from  the  feasible  set  of  parameter  values. 


6  =  max  log  (  dL  (8  j  z) )  (  14  ) 

N 
=  max  I   —  )   {)f     (i)B-i(i)  V(i)  +  log|B(i)j} 


[-71 


6f9       i=1 

(  15  ) 


The  log-likelihood  function  depends  on  the  innovations 
and  their  covariance.  To  optimize  the  likelihood  function,  a 
way  must  te  found  for  determining  these  quantities.  Both 
innovations  and  their  covariance  are  outputs  of  an  extended 
Kalman  filter. 
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The  extended  Kalman  filter  is  conventionally  divided 
into  two  parts.  In  the  first  part,  called  the  prediction 
equations,  the  state  equations  and  state  estimate  covariance 
equations  are  propagated  in  time  from  one  measurement  point 
to  the  next.  In  the  second  part,  called  the  measurement 
update  equations,  the  measurements  and  associated 
measurement  noise  covariances  are  used  to  improve  state 
estimates.  The  covariance  matrix  is  also  updated  at  this 
point  to  reflect  the  additional  information  obtained  from 
the  measurements. 

1 •   Prediction  Equations 

The  state  prediction  is  done  using  the  equations   of 

motion,  Appendix   B.  Starting   at   time   t      with   current 

i^1 

estimate  x (i-1 | i-1)  of  the  state  x (t  )   and   the   covariance 

i 

P(i-1|i-1),  the   following   equations   are   used  to  find  the 

predicted   state   x(i|i-1)   and   the   associated   covariance 
P(i|i-1)  • 


—   x(t|t.       )    =    f(x(t|t.       ),u(t),9,t)  (    16    ) 

dt  1-1  i-1 


T 
P(t|t.        )     =    F(t)     P{t|t.       )     +    P(t|t.       )     F     (t) 
1-1  1-1  1-1 


♦  rQfT  ( 1? } 


t    <  t  <  t 

i-1  i 
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The  n  x  n  matrix  F   is   obtained   by   linearizing   f 
about  the  best  current  estimate 


jf  (x(t|t.  )  ,u(t)  ,e,t) 

F(t)  = -7-77^ ; (  18  ) 


a*(t|t.   ) 
1-1 


Using  (  16  )-(  18  ),  we  can  obtain 


x(t  |t    )  £  x(i|i-1) 
i  i-1 


and 


P(t.  |t    )  A  P  (H  i-1)  (  19  ) 

l   1-1 


Thereafter,  the  measurement  update  equations  are  used. 

2.   Measurement  Update  Equations 

The  covariance  and  state  estimate  are  updated  using 
the  measurements.  The  necessary  relations  are  derived  in 
Eef.  4,  and  are  presented  here  without  proof.  The  innovation 
and  its  covariance  are 


V(i)  =  y(i)  -  h(x(i|i-1),u(t.),e,t.)        (  20  ) 

i     i 


T 
B  (i)  =  H  (i)  P(i|i-1)  H  (i)  +  R  (  21  ) 
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where  H  is  obtained  by  linearizing  h, 


dh(x(i|i-1)  ,u(t.,e,t.) 
H(i)  =  ■,.,.,.  ~ —  (  22  ) 

ax  (1 1 1-1) 


The  Kalman  gain  and  the  state  update  equations  are 


T 
K(i)  =  P(iji-1)  H  (i)  B-Mi)  (  23  ) 


x(ili)  =  x(i|i-1)  +  K(i)  V  (i)  (  2a  ) 


Finally,  F(iji),  the   covariance   of    error   in   updated 
state,  is  obtained  by 


P(i|i)  =  (I-K(i)  H(i))  P(iji-I)  (25) 


One   can   now   return   to   the   time   update   or   prediction 
equations. 

^ •   Optioization  Procedure 

Many  possible  numerical  procedures  can  be  used  for 
this  optimization  problem.  Modified  Newton-Raphson  or 
Quasilinear ization  have  been  found  by  experience  to  give 
quicker  convergence  than  most  procedures  such  as  the 
conjugate  gradient  or  the  Davison  method.  The  modified 
Newton-Raphson  method  is  a  second  order  gradient  procedure 
requiring  computation  of  first  and  second  order  partials  of 
the  log-likelihood  function. 
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a  log  ( se(eiz)) 

ae. 


N   r 
i=i  L 


VT(i)B-Mi)    *'M 


ae. 


1    /r  .    aB(i)  , 

-—1/       (l)B-Mi)- B-Mi)  V(i) 

2  a® . 


1        /  JB(i) 

2     I  ae 


D    /J 


(    26    ) 
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also 


32  log   (  *(etz)) 

3         * 


--I 


-B-»  (i) 


i-1 


a* 


ae 


T 

a  V    (i)  3B(i)  .       . 

b-1  (i)  -^ — B-i  (i)  V'(i) 


a  ^(D 

ae. 


a© 


.    aB(D         .     ,  . 

B-i(i)  — — B-»  (l)    V'(l) 

ae. 


ir/  T      .          ,     .     3B(i)        .     ..  aB(i)  .         .   . 

+  V  (i)  B-i  (l) B-i  (l) B-i  (l)    \/(i) 


ae 


ae. 


1  /  as(i)  aB(D 

Tr         3-i  (l)  B-i  (l) 

2  ae  ae 


a2^(i)  ... 

3      * 


•  T    •          .     •       d2B(i)  .        .    . 

-   V     (i)B-MD  -— — B-i  (i)  V/(i) 

ae.ae 


1       /  a2B(D 

+  — Tr         B-i  (l) 

2         I  ae .  ae 


j-    k/J 


(    27    ) 


3fK     —     \tZ,  ...,[) 


The  last  three  terms  in  the  equation  for  the  second 
partial  of  the  log- likelihood  function  involve  second 
partials  of  innovation  and  its  covariance.  These  terms  are 
usually  dropped        which        yields        the        second        partial 
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approximation 


32  log   (^e(eiz 
d«     ae^ 


N     P  T 

i-1  L  k 


LB-Mi)   ar(i) 


ae 


j 


8^   (i)  3B(i)  .      ,   . 

B-MD  "^ — B-i  <i)  v'(i) 


ae. 

3 


ae 


aB(D         .     .  . 

fl-Mi)  — ; b-i  (i)  V(i) 

ae. 


^  T      ■        ,    .    3B(i)       .    .    aB(i)      ,    .       , 

♦   V  (l)B-i  (1)— 3-1  (1)— B-i  (1)  V'(l) 


ae. 

3 


ae 


k 
3B(i) 


/         .    as(j) 

Tr        B-i(i)   B-»  (l) 

ae  ae 


k   /J 


(    28    ) 


From  equations  20  and  21  we  can  obtain  the  gradients 

of  innovation  and  its  covarience  for  parameter  0   as, 

j 


a  yd) 

ae. 


ah  |                jfx(iii-i) 
1x1    .,.,.,.      ae 

X  =  X  (1|  1-1) 

ah 
ae 

j    =    1,2,.. .  ,p 

i    =    1,2, ,N 

(  29  ) 


29 


3B<i)    3H(i)    .  .     T  .       .   dP(i|i-1)   T  . 
P(l|l-1)H  (1)  +  H(i)  ; H  (1) 


ae     39  ae 

3  3  3 

T 
.  .     P  (i)     3R 
+  H(i)  P(i|i-1)  ♦ (  30  ) 

ae.         ae. 

j      j 

j  =  1,2,. . . ,p 
i  =  1,2,. ..,N 


Recursive  equations  can  be  obtained  for  gradients  of 

the  predicted  state  and   its   covariance.  This   is   done   in 

stages   by   using   the   prediction   and   measurement   update 

equations  of  the   Kalman   filter.  Differentiating   equations 

16-18  with  respect  to  0   yields 

j 


ax(tjt.    )     af(x(tit.    ),u(t),t,e) 

d  i-1  1-1 


dt     a®  a® 


af(x(tit.    ),u(t),t,e)  ax(tit.    ) 
1-1  1-1 


+ 


ax(tit.    >  ae. 

1-1  3 


ax  (e) 
ax  (oio)       _o 

ae  a© 

j  j 


(    31    ) 
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3p(tjt.    ) 

d  i-1 

dt        39. 


-.t,,^  ap<tit.    ) 

3F(t  1-1 

— P(t|t.        )     ♦    F(t) 

ae                1-1  3e  # 


3P(t|t      )  : 

i-1       T  3F 

+  .  p     (t)     +    p(t|t.     J  — 


ae 


i-1     3e 


arQrT  +  r^rT  +  rQ3r 


ae. 


ae 


ae 


ae. 

3 


3P   (©) 
o 

ae. 


t        <  t  <t 

i-1 


(    32    ) 


g2f  (X(t|t       ),u(t),e,t) 

a^(t)  = i^i 

3e  39  3x(tjt.     ) 

j  j  1-1 


(    33    ) 


j    =    1,2,.  .  .  ,p 


The   sensitivity   functions,  used   in  the  SCIDNT  integration 

routines,  are    updated    at    measurement    points    by 

differentiating  equations  22  -  25  with  respect  to  9  . 

j 


i*h(z(i|i~1),ii(t.),e,ta) 

3H(1)  1      l 


ae. 

3 


39.3x(iji-i) 

3 


(  34  ) 
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3Mi) 


ae  f         a© 


aP(i|i-1)       T    . 

H     (1)     B-MD 


T 

.    .  3H    (i) 

+    P(lji-1) B-»{1) 

a© 


T                         .        3B(1) 
-    P(i|i-1)     H     (l)     B-i  (1)  — B-i  (i)      (    35    ) 

ae . 


ax  (Hi)      ax(iii-D      aK(D  ,  .         .a  v'a) 

-  +  —         -V(i)     +    K(i)— - (    36    ) 


ae 


ae 


ae 


3 


ae 


aP(ili) 

3e. 


=     (I-K(i)     H(i)) 


ae. 


aK(i) 
ae. 


H(i)    P(i|i-1) 


aH(D 

-    K(i)  -— P(iji-1) 

3e  . 

3 


(    37    ) 


j     =     1,2,3,... rp 


The  negative  of  the  matrix  of  second  partials  of  tae 
log-likelihood  function  is  called  the  information  matrix 
M.    The   step    size    A  6    for   parameter   estimates    is   given   by 


Ae  =m-* 


a  log    (^(eiz) 

ae 


(    38    ) 
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The   information   matrix  provides  a  lower  bound  on  parameter 
estimate  covariances,  i.e., 


T 
E  (9-9)   (9-6)   >  M-i  (  39  ) 


This  is  the  Cramer-Bao  lower  bound.  The   maximum   likelihood 
estimates  approach  this  bound  asymptotically. 

1  •   Lill§ar  Systems 

In   a   linear   system,  the   functions   f   and   h  are 
defined  as 


f(x,u,9,t)  £>F(9,t)x  +  G(9,t)u 


and        h(x,u,9,t)  £>H(6,t)x  +  D(9,t)u  (  40  ) 


For  the  linear  case,  it  has  been  shown  that  if  the 
model  is  correct,  the  innovations  are  white  and  nave 
Gaussian  density  at  the  true  values  of  the  parameters.  The 
assumption  of  fast  sampling  rate  is,  therefore,  not 
necessary, 

The  basic  algorithm  for  solution  of  the  linear  case, 

is   the   same,  however,  some   of   the   equations   can  now  be 

simplified.  The  equivalence  between  F(t)  and  H (i)  of   <  18  ) 

and    (  22  )    and    F(9,t)    and    H(9,t  )    of   (  40  )   is 

i 

obvious.  Equations  (  16  )  and  (  20  )  now  become 
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and 


—  x(t|t.   ) 
dt      i-1 


F(e,t)  x(t|t   )  +  G(e,t)  u(t) 
i-1 


(  41  ) 


V'(i)  =  Y(i)  -  H(9,tJ  x(i|i-1)  -  D(e,t.)  u(t  ) 

1  11 


(  42  ) 


Equations  (  29  )  and  (  31  )  can  be  written  as 


3H  (6,t  ) 

a/(i)  ax(i|i-1)  i   *  .  . 

=  -  H(6,t.)  ; ; x(i|i-1) 


ae 


i       a© 


ae. 

3 


|D(6,t.) 

; —  u  (t  . ) 

ae.  i 


(  43  ) 


and 


a     3x(t|t.  ) 

d        i-1 


dt   ae 


=  F(e,t) 


ax(tit.    ) 

1-1 

ae. 


jF(Qft)  A 

+ £ x(t|t.   ) 

ae  1-1 


3G(6ft) 

— pr-u(t) 

J 


(  44  ) 


All  other  equations  remain  the  same. 

There  is  a  considerable  reduction  in  the  computation 
requirement  for  time-invariant  linear  systems.  In  this 
case,  matrices  F,  G,  H,  D,  P ,  Q  and  R  and  their  derivatives 
with  respect  to  the  parameters  are  constant. 
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5«    Time  Invariant  Linear  Systems  In  Statistical  Steady 
State 

In  many  aircraft  applications,  the  Kalman  filter  is 
in  steady  state  for  the  duration  of  the  experiment.  This 
occurs  when  the  Kalman  filter  is  in  operation  for  a 
sufficiently  long  time  and  the  process  and  measurement  noise 
covariances  are  constant.  The  Kalman  gain  and  the 
innovations  and  the  state  covariances  approach  constant 
values.  The  time  update  and  measurement  update  equations  for 
the  covariances  are 


—-p(tit.  )  =  FP(tit.  )  +  p(tit.  )f  +  for 

dt      1-1  1-1  1-1       '    ■ 


By  definition  of  the  steady  state 


(  45  ) 


T 
K  =  P(i|  i-1)H  B-i  (  46  ) 


T 
B  =  HP  (i|i-1)  H   +  R  (  47  ) 


P(i|i)  =  (I-KH)  P(i|i-1)  (  48  ) 


P(i-1ji-1)  =  P(i|i) 


Therefore,  from  (  45  ) 
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T 
F    At  F       At 

P  (i|i-1)    =    e  P(i-1|i-1)    e 


t 


e        1  Q|      e  1  dT 


i-1 


T 
F    At  F      At 

=    e  (I-KH)     P(i|i-1)e 


/F(t.-T 


-T)      _     _T    F     (t    -T) 

PqP    e  1         'dT 


i-1 


A^(At)      (I-KH)     P(i|i-1)    0(At)     +    Q« 

T 
=    j*(At)      (P(i|i-1)     -    KBK    )     0(At)     +    Q» 

(    49    ) 


Using    (    46    )  ,     (    47    )    and    (    49    ),    we   can    solve      for      P(i|i-1) 
and    then    find   K    and    B.    Also,    it   can   be    shown    that 


5P  3P         T  T 

=    A     A      +    A      -   #PA    P0 

36  1     ,)e  1  2  3 


(    50    ) 


i)K                              AP         T  T 
=     (I-KH)  H        (HPH       +    R)"i     +    A 

ae .  a© .  2 


(    51    ) 
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iB 

an      t        ap     t         an       aR 

PH   +  H  H   «•  HP  + 

ae. 

ae               a©              a©      ae 

3 

3           3           3      3 

(  52  ) 


where 


A   =  0(I-KH) 

1 


30  T  10 

A   =  (I-KH)  20       +  0(I-KH)P 

2   ae 


ae 


3H    t   jq« 

-  0K  -T—  P0   + 


ae  ae 

3         3 


3H  T     3H     T 

A   =  B"lH  4-  H  B-1  PH 


3   ae 


ae 


aH '     a* 

+  HP  ♦  B-*H 

ae  .     ae  . 

3      3 


(  53  ) 


P&P(i|i-1) 


ap 

Thus,  it  is  possible  to  solve  for  using  (  50  )  and  then 

ae. 

find  and  -  from       (    51    )       and       (    52    ) .    Equation       (   50    ) 


ae  ae 

3        3 
is   a   linear   equation   in 


11 

ae 


and   the   coefficient  of  the 


unknown   matrix    does   not    depend   on    the    parameter 
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6  .  Thus,  the  sensitivity  of  the  state  covariance  matrix  can 

J 

be  determined  very   quickly   for  all   parameters.  Once  the 

sensitivity   of   P,  K  and   B   for   unknown   parameters  is 

determined,  only  state   sensitivity   equations  need   to  be 

updated.  The  computation  of  state  sensitivity  functions  can 

be  reduced  to  many  fewer  equations. 

An  approximation  simplifies  the  problem  further.  The 

unknown  parameters  are  defined  to  include  elements  in  K  and 

B  matrices  instead  of  Q  and  R.  Optimizing  the  log-likelihood 

function  for  parameters  in  B  gives 


1  t**^  T 

B  =  —  )      V'(i)  /  (i>  {  54  ) 

i=1 


The  gradient  of  the  log-likelihood  function  with  respect   to 
other  unknown  parameters  is 


j  1=1  j 


The   sensitivity   of  innovations  to  parameters  is  determined 
using  the  following  recursive  equations 


—  x(t|t.   )  =  Fx(t|t.   )  +  Gu(t) 
dt      1-1  i-1 
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ax(tit.    )  axttit.    ) 

d                     1-1            9F      *                                                 1-1 
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dt     ae  a®  i-i  a® 
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ii 
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u(t) 


j    =    1,2, ...,p 


t  <    t   <    t 

i-1  i 


(   56   ) 


V'(i)    =    y(i)    -   Hx(i|i-1)    -    Du(t.) 

i 


x  (iji)    =    x  (i|i-1)     +    K  v'(i) 


(    57    ) 


a  v'(i) 

2« 


an    .,.,.„      ,  ax(iii-i)      ao 
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ae.  ae  ae         i 


(    58    ) 


ax  (iji) 


ax(iji-i)      aK     ,  av'a) 

-  +  ■  v^(i)    +    K  — 


ae 


ae 


a« 


j    =1,2, 


Note   that 


=0  if  9      is   not   an    element   of   K    matrix 

de.  j 


(    59    ) 
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where  I     is  a  matrix  of  all  zeroes   except  a   1   at  the 

j'k« 
j1,  k*  position. 

This   approximation   simplifies   the   optimization 

considerably.  However,  this     usually     leads    to    an 

overparameterized  model.  In  other  words,  once  K   and   B  are 

determined,  it   is   not   possible   to  find  any  corresponding 

P,  Q   and   S   which   have   the   desired   (known   a   priori) 

structure.  A   good   estimate   of   elements   in   P ,  Q   and  R 

matrices   can  be   obtained   by   a    least-squares    type 

approach.  The  fit   to  the  observed  data  is  better  than  with 

true  values  of  parameters  but  the  parameter  estimates  do  not 

have   minimum   variance.  This   approximation   is  not  good  in 

aircraft  application  where  the  structure  of  Q  and  R  is  known 

fairly   well,  but   is  excellent  where  there  are  many  process 

noise  sources  and  the  characteristics  of  both  Q   and   R  are 

relatively  unknown  (e.g.,  economic  systems). 

6«   U^xiEum  Likelihood  With  No   Process   Or   Measurement 
Noise 

The  maximum  likelihood  method  can  be  simplified  when 
either  process  noise  or  measurement  noise  are 
absent.  Assuming  no  process  noise  we  can  proceed  as 
follows. 

If  the  process  noise  is  zero  and  initial  states  are 
known  perfectly,  i.e.,  w(t)  and  P(0)  are  zero,  the 
covariance  of  the  error  in  the  predicted  state  is  also 
zero.  It  is  clear  from  (  23  )  that  the  Kaiman  gains  are 
zero.  The  innovations  are  the  output  error,  i.e., 


\T(i)  =  y(i)  -  h(x(t.),u(t.),6,t.)  (  60  ) 

111 
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and  the  innovation  covariance  is  (  21  ) 


B(i)  =  R 


(  61  ) 


the  log- likelihood  function  is, 


log  (se(6|z))  = 


n 

1  T" 

~yH  ^  (i)  R_iv/(i)  +  iog  ,R| 


i-1 


(  62  ) 


which  on  optimizing  for  unknown  parameters  in  R  yields 


R 


=  —  y  V'(i)   V  (i) 
i=1 


(  63  ) 


The  equality  in  (  63  )  holds  only  for  those  elements  of  R 
which  are  not  known  a  priori.  For  instance,  even  if  R  is 
known  to  be  diagonal,  the  right  hand  side  matrix  will  not  be 
diagonal  in  general,  but  the  off-diagonal  terms  should  be 
ignored  before  they  are  equated  to  R.  Using  (  63  )  in  (  62  ) 
we  obtain 


log  (c*(0|z))  ="—  ]T  ^  (i)  R_1  V(i) 


i=1 


♦  constant 


(  64  ) 
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The  optimizing  function  is  the  same  as  that  for  the  output 
error  method  except  that  the  measurement  noise  covariance 
matrix  is  determined  using  (  63  )  and  is  used  as  the 
weighting  matrix  in  the  criterion  function.  In  the  output 
error  method,  the  measurement  noise  is  assumed  known  and  the 
weighting  function  is  arbitrary. 

The  first  and  second  derivatives  of  the 
log-likelihood  function  with  respect  to  unknown  parameters 
are 


ae 


log(3P(6|z)J    = 


i=1 


V     (i)     R-i  -i— - (   65    ) 
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(    66    ) 


The  terms  in  the  second  derivative  are  approximated  as 


a?log(  £{Q\z)) 

&  .   aev 
1        k 
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i=1 
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-  R-i 


ae, 


ae. 
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(  67  ) 


With  no  measurement  noise  we  can  proceed  as 
follows.  If  all  states  are  measured  with  no  noise,  the 
covariance   of   the   error  in  state  estimates  is  zero  at  the 
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beginning  of  any  time  update, 


P  (i-1|i-1)  =  0 

(  68  ) 


and       x(i-1|i-1)  =  x(i-1) 


It  is  easy  to  show  in  this  case  that  for  fast  sampling  the 
log-livelihood  function  is  quadratic  in  the  difference 
between  measure  values  of  x  and  f(x,u,6,t).  The  method 
reduces  to  the  equation  error  method,  the  weight  W  being 
chosen  as 


1     r  .  .  t 

W  = —  /  (x  -  f(x,u,6,t))  (x  -  f  (x,u,9,t))  dt 
T  J 


(  69  ) 


Thus,  the  maximum  likelihood  method  and  equation-error 
methods  are  equivalent  except  for  the  technique  for  choosing 
the  weighting  matrix. 


F.   AIRCRAFT  PARAMETER  IDENTIFICATION  USING  THE  MAXIMUM 
LIKELIHOOD  METHOD 

In  applying  the  maximum  likelihood  theory  it  is 
necessary  tc  consider  the  aircraft  equations  of  motion  in 
discrete  form.  It  was  assumed  that  the  structure  of  the 
model  was  known.  The  vector  of  unknown  parameters  from  x  (0)  , 
#*      G,  P  ,   H   ,   Q,   and    R    was   denoted   by    6.  The 
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problem,  therefore,  was  to  identify  9  using  the  maximum 
likelihood  method.  In  constructing  the  problem  it  was 
further  assumed  that  the  aircraft  which  had  been  flying  for 
seme  time  was  in  a  steady  state  condition,  and  process  noise 
(wind  gust  effects)  was  not  present. 

In  the  computer  implementation,  the  negative 
log-likelihood  function  was  minimized  using  a  Gauss-Newton 
gradient  procedure.  A  least  squares  routine  was  utilized  to 
provide  approximate  airframe  parameter  start-up 
values.  Additionally  the  results  of  the  SCIDNT  computer 
program  were  utilized  to  obtain  frequency  and  damping. 


G.   SCOPE  OF  RESEARCH 

Research  was  limited  to  the  longitudinal  case  since  both 
available  computer  and  research  time  were  limited.  Although 
additional  data  were  desired,  the  project  airplane  was 
substantially  damaged  by  fire  and  was  removed  from  further 
test  activity.  Further  data  gathering  with  the  F-14A 
airplane  was,  therefore,  not  possible.  Specific 
longitudinal  stability  derivatives  identified  are  given  in 
the  Results  and  Discussion  section.  There  were  no  airplane 
flight  restrictions  which  affected  the  test  methods  or 
limited  the  scope  of  research.  The  research  was  limited  to 
establishing  base  line  data  for  the  short  period  and  phugoid 
modes  and  thence  to  determine  the  feasibility  of  identifying 
these  modes  simultaneously. 


H.   METHOD  OF  TESTS 

In  general  test  proceedures  were  in  accordance  with 
those  contained  in  the  U.  S.  Naval  Test  Pilot  School  Fixed 
Wing  Stability  and  Control  Manual,  Ref  5.  There 
were,  however,  some  non-standard  test  techniques  utilized  in 
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exciting  the  desired  airplane  modes.  The  short  period  mode 
was  excited  by  the  standard  doublet  method  and  additionally 
by  applying  a  sinusoidal  control  input  at  varying 
frequencies.  The  phugoid  mode  was  excited  by  the  Delta 
Airspeed  method  in  all  cases,  [5]. 

In  exciting  the  phugoid  and  short  period  modes 
simultaneoulsy,  a  combination  of  standard  test  techniques 
were  utilized.  The  "simultaneous"  technique  consisted  of 
establishing  the  airplane  in  a  phugoid  maneuver  and  during 
the  first  quarter  cycle  of  the  maneuver  applying  a  doublet 
input  to  excite  the  short  period  mode. 

The  data  system  utilized  within  the  airplane  was 
specifically  designed  to  be  compatible  with  the  Naval  Air 
Test  Center  Real-Time  Telemetry  Processing  System  (RTFS) .  A 
brief  description  of  the  airplane  and  ground  data  systems  is 
presented  in  Appendix  C.  Data  reduction  was  performed  en  the 
RTPS  Xerox  Sigma  9  computer  in  the  batch  processing 
mode.  Raw  airplane  data  were  converted  to  an  Engineering 
Unit  tape  and  used  as  input  to  the  digital  parameter 
identification  program  (SCIDNT) .  In  addition  a  least  squares 
parameter  identification  routine  was  utilized  to  generate 
parameter  start  up  values  for  SCIDNT.  This  least  squares 
routine  extracted  a  set  of  approximate  stability  derivatives 
by  minimizing  the  error  between  the  model  and  actual 
airplane  state  vector  time  histories. 
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II.   RESULTS  AND  DISCUSSION 


A.   SHOEI  PERIOD  AIRFRAME  PARAMETER  IDENTIFICATION 

1 •   Air  plane  T est  Conditions 

Initial  trimmed  test  conditions  for  each  control 
input  are  presented  in  Table  I.  Prior  to  each  control  input 
the  airplane  was  trimmed  for  steady  state  "hands  off"  flight 
[5]. 


Altitude 
(ft  HPC) 

MACH 

CG 
(%  MAC) 

Loading 

29  925 

.6149 

11.18 

C  1) 

B 

30  503 

.6289 

10.985 

B 

Sine  Wave  Input 
Doublet  Input 


(1) .   Four  Sparrow  III  missiles  (AIM-7)  on  stations  3,  4,   5 
and  6 

lable  I.   Short  Period  Test  Conditions 


2.   Method  Of  Investigation 

The   airplane   model  equations  used  were  those  shown 

in  Appendix  B.   In  order  to  isolate  the   short   period   mode 

from   the   longitudinal  equations  of  motion  the  state  vector 

parameters  oc  ,  q,  and  8  were  made  insensitive  to  u.  Thi^s  was 

accomplished   simply   by   fixing  X   ,  x  ,  X  ,  Z  ,  M   and  X^ 

u    q    u    u       e 
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equal  to  zero  and  not  identifying  any  of  thera.  This  equation 

modification   resulted    in    the   clasical  short   period 
approximation  Ref.  6. 

To  provide  SCIDNT  with  approximate  airframe 
parameter  start-up  values, the  least  squares  routine  was 
utilized.  For  the  short  period  case  the  least  squares  start 
up  values  were  very  good,  a  comparison  of  a  typical  set  of 
least  squares  parameters  and  the  final  SCIDNT  parameter 
values  is  shown  in  Table  II. 


Z  Z        MM  Z~  M-, 

Least       -.4909  0.9624  -1.78    -.9085  0.210  6.09 
Squares 

SCIDNT       -.4903  0.9805  -1.95    -.795  0.278  5.542 


Table  II.   Comparison  of  Least  Squares  Parameter 
Values  and  SCIDHT  Parameter  Values 

In  actual  practice,  however,  airplane  data  was 
recorded  in  long  bursts  which  usually  included  three  or  more 
maneuvers  ,  executed  with  near  identical  initial 
conditions.  As  a  result  once  a  set  of  start  up  values  were 
determined  for  a  given  test  condition  it  was  unnessary  to 
determine  a  set  of  start  up  values  for  other  like 
maneuvers. 

The  actual  parameter  identification  proceedure  ,  for 
the  short  period  mode,  was  relatively  simple  once 
familiarity  with  the  digital  program  was  attained.  Upon 
obtaining  a  set  of  "reasonable"  initial  parameter 
estimates,  all  six  short  period  parameters  were  estimated  in 
an  iterative  fashion.  Here  "reasonable"  denotes  proper 
parameter  sign  and  magnitude.  Each  computer  run  was  made 
with  the  number  of  iterations  and  step  cuts  both  initialized 
at  ten.  During  a  given  iteration  the   program   attempted   to 
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minimize  the  likelihood  function.  If  the  likelihood  function 
increased/  a  series  of  step  cuts  would  be  taken  until  the 
likelihood  function  again  decreased  or  the  maximum  number  of 
step  cuts  were  exceeded.  The  program  would  then  continue  to 
iterate  until  the  desired  number  of  iterations/step  cuts 
were  reached.  Due  to  program  logic  one  more  step  cut  than 
requested  was  executed.  For  the  short  period 
mode,  however,  convergence  was  obtained  early  in  program 
execution.  In  general,  iterations  conducted  with  doublet 
control  input  data  converged  very  quickly  and  began 
executing  step  cuts  until  the  allowable  number  was 
exceeded.  Figure  2  shows  the  change  in  the  likelihood 
function  during  iteration  for  a  doublet  input.  Program 
convergence  had  been  essentially  reached  by  the  fourth 
iteration  at  which  time  a  series  of  step  cuts  were  made. 
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Figure  2.   Likelihood  Function  Change  During  Iteration  For 

Doublet  Input 


The   convergence   characteristics   using   sinusoidal 
input   data   were   somewhat   different.  During   the  program 
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iteration  the  likelihood  function  continuously  decreased  and 
no  step  cuts  were  taken,  Figure  3.  Program  convergence  for 
this  case  was  evidenced  by  a  leveling  out  of  the  likelihood 
function,  but  not  until  the  eighth  iteration.  It  was  felt 
that  since  the  data  length  for  the  sinusoidal  input  was 
longer  than  the  short  period  data  length  the  larger  number 
of  iterations  until  program  convergence  was  not  abnormal. 
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Figure  3.   Likelihood  Function  Change  During  Iteration  For 

Sinusoidal  Input 


After   convergence   was   obtained   for  the  parameter 

estimates,  the   Angle   Of   Attack   Scale   Factor   (K   }   was 

determined.  K  ,  was   a   scaling  factor  in  the  angle  of  attack 
cC 

measurement   equation,  Appendix   B.  In   estimating   K     two 

terms   in   the   model  were  actually  affected,  K^  itself  and 

the  variable  (1   K   /U)  which  was  used  to  multiply  q  in   the 

measurement   equation.  Since,  in  essence,  two   terras   were 

being  estimated  a  modified  form  of  parameter  estimation   was 

utilized.  In   the   term   (1_  K   /u")   the   value   of   K    was 

*  cC  t£   '      '  CC 


initially   set   to   1.0   and   the   single    term 


K 


oc 


was 
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determined.  This   new  value  of  K  ,  was  then  substituted  into 

the  term  (1   K  /U)  and   again   a   set   of   iterations   were 

performed   to   determine   an   updated   K   .  This   method   of 

iteration  was  continued   until   there   was   essentially   no 

change   in   K   ,  at   which   time   the  value  was  said  to  have 

converged.  In  most  cases  only  two  iterations   were   required 

to  obtain  a  converged  K^.  . 

Upon  determining  a  value   for   K   ,a   final   set   of 

iterations    were  •  performed    to   determine   values   for 

Z,X,M,6.At  this  point,  the  short   period   mode   time 

o   o    o    o 
history    fits   were   nearly   identical   and   the   estimate 

statistics  indicated  a  high  degree   of  confidence.  Attempts 

at   further   iterations  to  improve  the  estimates  resulted  in 

using   a   lot   more   computer   time  ,  with   very    little 

improvement.  Therefore,  the   procedure  utilized  to  estimate 

the  short  period  mode  was  to  first  estimate  the  short  period 

parameters,  second   estimate   K   ,  and  finally,  estimate  the 

initial  values. 

3 .   Data  Analysis 

Data  from  typical  short  period  estimation  runs  for  a 
doublet  and  sinusoidal  input  are  presented  in  Tables  III  and 
IV.  Typical  short  period  time  history  fits  are  presented  in 
Appendix  D  Figures  1  and  2,  for  the  doublet  and  sinusoidal 
inputs  respectfully.  It  should  be  noted  that  the  scales  for 
the  short  period  time  histories  are  different  for  the 
doublet  and  sinusoidal  inputs.  Furthermore,  since  the 
estimated  stability  derivatives  were  in  dimensional  form  the 
actual  derivative  values  differed  slightly  from  maneuver  to 
maneuver . 

The  data  obtained  with  the  doublet  input  in  most  all 
cases  converged  to  believable  derivative  estimates,  i.e., 
expected  sign  and  magnitude,  with  favorable  iteration 
statistics.  For  example,  refering  to  Table  III,  the  value  of 
M_  was  -1.84*1  with  a  standard  deviation  of  0.0061  and  a 
final   step   size  of  -0.0146  which  is  significantly  lower  in 
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magnitude  than  the  estimated  parameter.  From  probability 
theory  we  can  say  that  the  probability  is  99.74%  that  the 
value  of  M^  is  -1.844  ±0.0183  which  is  well  within 
acceptable  engineering  tolerances.  This  same  analysis  was 
carried  out  for  the  other  estimated  derivatives  with  similar 
results. 


Derivatives 

Z 

OC 

Z 

g 

M 

OC 

M 

g 

z*e 

M*e 

Estimate 

-.4138 

1.099 

-1.844 

-, 

.826 

0.1628 

5.098 

Standard 
Deviation 

0.0024 

0.0038 

0.0061 

0, 

.0037 

0.0055 

0.0104 

Parameter 
Step  Size 

-.0006 

-.012 

-.0146 

-  , 

.0001 

0.0143 

0.028 

Eigenvalue 
Spread 

1.28 

X 

10& 

Likelihood 

Function 

Value 

-5.  1 

X 

103 

Time 

History 

Matches 

CT   =  0, 

oC 

.00396 

cr    = 

u 

0. 

.564 

cr  =   o.i 
q 

3039 

cr   =  o.oo43  cr  =  2.33 
e  n 

z 

Table  III.   SCIDNT  Estimates  of  Longitudinal  Short  Period 
Stability  Derivatives  Using  a  Doublet  Input 


The  eigenvalue  spread,  that  is,  the  difference 
between  the  largest  eigenvalue  and  the  smallest  eigenvalue 
of  the  information  matrix,  was  1.28  x  106  which  was  less 
than  the  maximum  spread  allowed  by  SCIDNT  (1010).  All 
parameters  converged  indicating  no  identif iability 
problems,  and  the  eigenvalue  spread  was  therefore  not 
considered  significant.  The  likelihood  function  value  was  a 
large  negative  number  which  converged  rapidly  (Figure 
2) ,  further  indicating  no  identif iability  problems. 
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As  a  final  check  of  stability  derivative  correctness 
a  series  of  time  history  fits  were  constructed.  In 
constructing  these  plots  both  the  model  and  airplane  were 
excited  with  the  same  control  input  and  the  airplane  time 
history  plots  superimposed  on  the  same  scale.  The  time 
history  fit  for  the  doublet  input  was,  in  general,  very 
good,  with  the  difference  between  the  measured  value  and 
estimated  value  within  the  instrumentation  measurement  error 
tolerances  Ref.  7. 

Referring  to  Table  IV,  it  can  be  seen  that  typically 
the  estimated  derivatives  using  a  sinusoidal  input  converged 
tc  appropriate  (sign  and  magnitude)  stability  derivative 
values. 


Derivatives 

Z 

Z 

g 

oc 

M 

q 

Z<f 
e 

l\ 

Estimate 

-.4303 

1.  1229 

-2.007 

-.7615 

0.1532 

5.2274 

Standard 
Deviation 

0.0042 

0.0061 

0.013 

0.007 

0.0045 

0.0181 

Parameter 
Step  Size 

0.0007 

-.0016 

-.0043 

-.0061 

-.0037 

0.0332 

Eigenvalue 

1.67 

X   105 

Spread 

Likelihood 

Function  -7.323  x103 

Value 

Time 

History  CT  =    0.00533  CT     =    0.9  <T     =   0.006 

Matches         *  u  q 

CT     =  0.01    0~    =    2.32 
e  n 

z 

Table  IV.   SCIDNT  Estimates  of  Longitudinal  Short  Period 
Stability  Deravitives  Using  A  Sinusoidal  Input 


In  general,  however,  the  parameter  standard  deviations  and 
step  sizes  were  not  as  favorable  as  those  obtained  with  the 
doublet     input     for     these     two    examples.  It    was 
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possible,  however,  to  select  another  pair  of   maneuvers   and 

reach   the   opposite   conclusion,  ie.  sinusoidal   input  more 

favorable.  Since   all   parameter   estimates   for   both    the 

doublet   and   sinusoidal   inputs   converged  it  was  concluded 

that  the  slight   difference   in   statistical   data   was   not 

significant.  If  a  determination  as  to  which  control  input  to 

use   were    required,  it    was    felt    that    the    primary 

consideration   would  be  ease  of  pilot  input.  If  an  autopilot 

device  were  available  which  could  input   the   sinusoid   then 

this  input  should  be  used,  however,  with  no  autopilot  system 

installed   the   doublet   input    would    provide    excellent 

results. 

A  consistent  problem   that   occured   throughout   the 

research   was  concerned  with  the  initial  conditions  prior  to 

the  control  input.  If  the  initial  condition  of  the   airplane 

was   not   stabilized  prior  to  the  control  input  or  if  a  data 

burst  was  taken  in  the  middle  of  a   control   input,  no   time 

history   match   was   possible.     According   to  theory  it  is 

possible  to  extract  the   airframe   parameters  at   any   time 

during   a   maneuver,  this,  however  did   not   appear   to   be 

true.  Figure  3  Appendix  D  shows  the   results  of   estimating 

the     parameters    starting    from     non-zero     initial 

conditions.  This  time  history  fit  was   generated   after   all 

iterations   on   both   parameters  and  initial  conditions  were 

completed.  The  parameter  estimates  did  not  converge  and   the 

most   powerful  derivative,  H   ,  was  much  larger  in  magnitude 

cc 

than  any  ether  estimates  with  a  value  of  -3.7.  These  results 
occured  any  time  an  attempt  was  made  to  start  SCIDNT  from  an 
initial  condition  other  than  steady  state.  It  was  concluded 
that,  within  the  scope  of  this  research,  meaningful 
parameter  estimates  were  only  obtainable  when  started  from 
an  initial  steady  state  condition. 

A  comparison  of  short  period  characteristics 
obtained  from  Navy  Preliminary  Evaluation  (NPE)  data  and 
those  determined  by  SCIDNT  is  presented  in  Table  V. 
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Description 


Damping 
Ra 


pm< 
tio 


Frequency 


NPE 

SCIENT-Doublet 

SCIDNT-Sinusoidal 


C 


sp 


0.35 

0.403 

0.371 


n 
sp 
(rad/sec) 

1.49 

1.54 

1.61 


Table  V.   Comparison  of  Longitudinal 
Short  Period  Characteristics 


It  was  concluded  that,  within  the  scope  of  this 
research,  the  parameter  identification  program  SCIDNT  was 
easy  and  straightforward  to  use  in  determining  the  short 
period  mode.  In  addition,  it  was  concluded  that  it  was  not 
possible  to  obtain  accurate  estimates  when  starting  from 
non-steady  state  initial  conditions. 


B.   PH0GOID  AIRFRAME  PARAMETER  IDENTIFICATION 

1 •   Airplane  Test  Conditions 

Initial  trimmed  test  conditions  for  the  phugoid 
maneuver  are  presented  in  Table  VI.  Prior  to  the  control 
input  the  airplane  was  trimmed  for  steady  state  "hands  off" 
flight  [5], 


Phugoid  Input 


Altitude 
(Ft  HPC) 

30  055 


MACH 


0.614 


c  c 

(%    MAC) 
10.94 


Loading 
B 


Table  VI.   Phugoid  Test  Conditions 
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2.   Method  Of  Investigation 

The  airplane  model  equations  used  were  those  shown 
in   Appendix   B.  The   derivitives   X   and  X^  were  fixed  at 

zero   values   since   their   effect   on   the   phugoid   was 

small.  Additionally,  it   was   felt,  fixing  these  terms  would 

facilitate   the   identification   of   the   dominant    phugoid 

terms.  Finally   from  Appendix  B  it  can  be  seen  that  the  terra 

X7  was  defined  as  {X   -  oc  U  )   and   therefore   due   to   its 
q  q      O  o 

relatively    large    size   was   approximated   by   oC  0   and 

o  o 

fixed.  The  parameters  which  were  identified  during  the   pure 

phugoid  maneuver  were  therfore,  Z  ,  X  ,  and  M  . 

u   u       u 

Since  there  was  no  short  period  input  during  the 
pure  phugoid  maneuver  it  was  not  possible  to  effectively 
identify  the  short  period  parameters.  Additionally  the 
program  solution  diverged  without  a  set  of  predefined  short 
period  parameters  and  it  was  necessary  to  obtain  parameters 
from  a  prior  run.  A  set  of  short  period  parameters 
were,  therefore,  obtained  from  a  previous  identification  run 
and  fixed  throughout  the  phugoid  parameter  identification 
runs. 

In  an  attempt  to  provide  SCIDNT  with  approximate 
airframe  parameter  start-up  values  the  least  squares  routine 
was  utilized.  For  the  phugoid  case  the  least  squares  start 
up  values  were  unusable.  For  example,  least  squares  values 
for   X     was   -2.6697  and  X^  was  -240.5,  both  of  which  are 

terras  of  small  absolute  value.  As  an   additional   check,  the 

least  squares  parameter  estimates  were  used  to  generate  a 
set  of  time  histories.  For  these  time  historys  the  standard 
deviation  of  airspeed  and  normal  acceleration  was  49.12 
ft/sec  and  5.44  ft/sec2  respectively.  Therefore,  since  the 
least   sguares   routine   did  not  provide  acceptable  start-up 
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values  the  identified  parameters  were  set  initally  to  zero 
for  the  first  SCIDNT  iteration.  Further  attempts  to  obtain 
phugoid  start  up  values  for  SCIDNT  from  the  least  squares 
routine  proved  futile  and  its  use  was  discontinued  for  the 
phugoid  case. 

The  version  of  SCIDNT  utilized  for  this 
investigation  would  accept  a  maximum  of  500  data  samples  per 
run.  Since  the  phugoid  period  was  approximately  100  sec  a 
data  burst  recorded  at  20  samples  per  second  (SPS)  would 
only  include  the  first  25  sec  of  the  maneuver.  This  time 
span  was  obviously  acceptable  for  the  short  period  mode  but 
clearly  unacceptable  during  a  phugoid  maneuver.  An  attempt 
was  therefore,  made  to  utilize  data  at  2  and  5  SPS  to 
include  as  much  of  the  phugoid  maneuver  as  possible.  This 
data  was  from  the  same  maneuver  with  data  tapes  recorded  at 
2  and  5  SPS. 

Although  the  2  SPS  data  would  provide  a  longer 
phugoid  time  span  in  theory,  the  program  exhibited  severe 
instabilities  during  the  control  input  and  therefore  only 
the  latter  portion  of  the  data  burst  was  usable.  It  was 
strongly  suspected  that  these  instabilities  were  due  to  the 
unusually  rough  control  inputs. 

The  actual  parameter  identification  procedure  for 
the  phugoid  mode  was  somewhat  different  than  that  utilized 
during  the  short  period  investigation.  The  short  period 
stability  derivatives  were  obtained  from  a  run  whose  initial 
conditions  most  closely  approximated  those  of  the  phugoid 
run.  These  short  period  parameters  were  then  fixed  as 
constants  during  the  phugoid  iterations.  Computer  runs  were 
made  with  the  number  of  iterations  and  step  cuts  both 
initialized  at  ten.  The  program  would  then  continue  to 
iterate  until  the  desired  number  of  iterations/step  cuts 
were  completed.  For  both  the  5  SPS  and  2  SPS  data 
convergence  was  obtained  prior  to  executing  the  full  10 
iterations.  Iterations  conducted  with  phugoid  data  recorded 
at   2   SPS   appeared   to   converge   very   quickly   and  began 
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executing  step  cuts  until  the  allowable  number  was 
exceeded.  As  in  the  short  period  case  the  program  executed 
one  step  cut  more  than  requested.  Figure  4.  presents  the 
change  in  the  likelihood  function  during  iteration  using 
phugoid  data  recorded  at  2  SPS.  The  program  converged  to  its 
"best"  solution  by  the  third  iteration  at  which  time  the 
step  cuts  were  executed. 
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Figure  4.   Likelihood  Function  During  Iteration 

Using  2  SPS 


Note  that  although  the  program  had  apparently 
converged  at  the  third  iteration  the  value  of  the  likelihood 
function  decreased  until  the  fourth  step  cut.  The  function 
value  then  remained  essentially  constant  throughout  the 
remaining  step  cuts. 

The  convergence  characteristics  using  data  recorded 
at  5  SPS  were  somewhat  different.  Figure  5. 
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Figure  5.   Likelihood  Function  During  Iteration 

Using  5  SPS 


It  should  be  noted,  however,  that  the  scales  are 
different  for  each  plot.  The  program  had  essentially 
converged  at  the  fifth  iteration,  at  which  time  the  program 
executed  step  cuts.  Note  that  the  likelihood  function 
values  decrease  using  5  SPS  data,  and  were  essentially  the 
same  using  the  2  SPS  data.  However,  the  time  history  matches 
for  each  set  of  data  were  significantly  different.  The  angle 
of  attack  scale  factor  (K  )  was  not  estimated  since  the 
value  of  K  from  the  short  period  data  was  utilized.  After 
convergence  was  obtained  for  the  parameter  estimates  a  final 
set  of  iterations   was   performed   to   estimate   values   for 


Z  ,  X  ,    M  ,  6     and 
o    o    c    o 


the     three    phugoid    parameters 
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simultaneously.  At  this  point,  the  phugoid  time  history  fits 

were  acceptable  for  the  5  SPS  data  and  clearly  unacceptable 
using  the  2  SPS  data.  Further  attempts  to  improve  the 
estimates  and  time  histories  resulted  in  using  a  lot  more 
computer  time  with  no  improvement. 

3.   Data  Analysis 

Data  from  phugoid  estimation  runs  using  2  SPS  and  5 
SPS  are  presented  in  Tables  VII  and  VIII 
respectfully.  Phugoid  time  history  fits  are  presented  in 
Appendix  D  Figures  4  and  5  for  2  and  5  SPS  respectfully.  It 
should  be  noted  that  the  scales  of  the  phugoid  time 
histories  for  2  and  5  SPS  are  different. 

The  data  obtained  using  a  sample  rate  of  2   SPS   was 

totally    unacceptable.  For    example,  refering    to   Table 

VII,  the  value  of  H   was  -0.00000017   while   the   last   step 

u 

size  was  0.003. 


Derivatives     Z       X        X       X       M       X^ 

u       OC       u       g       u       °e 

Estimate     -.003      0      -.758   -47.9   -1.7x10-?   0 

Standard      .00006   fixed    .0542    fixed    .0002     fixed 
deviation 

Parameter    -.001    fixed    0.444    fixed   0.003      fixed 
Step  Size 

Eigenvalue  3.54  x  108 

Spread 

Likelihood 

Function  ^8.4  x  103 

Value 

History  0~  =    0.0048  0"      =   26.77      <T    =    0.00534 

Matches  OC  u  q 

0~    =    0.053      <T    =    2.53 
6  n 

z 

Table  VII.   SCIDKT  Estimates  of  Longitudinal  Phugoid 
Stability  Derivatives  Using  2  SPS  Data 


Although  the  eigenvalue  spread  and  likelihood 
function  change  during  iteration  indicate  no  significant 
program  difficulties  the  standard  deviations  of  the  state 
error  vector  are  very  large.  As  an  approximation  it  can  be 
said  that  the  phugoid  maneuver  is  essentially  a  varying 
airspeed  and  pitch  attitude  maneuver  with  a  constant  load 
factor.  From  Table  VII  it  can  be  seen  that  the  standard 
deviation  for  airspeed  was  26.77  ft/sec-1  and  for  pitch 
attitude  was  0.053  rad,  both  of  which  were  unacceptable. 

The  time  history  plots  constructed  using  the  2  SPS 
data  showed  no  real  agreement  between  the  model  and  the 
airplane.  In  constructing  these  plots  both  the  model  and 
airplane  were  excited  with  the  same  control  input  and  the 
airplane   time   history   plots   superimposed   on   the   same 
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scale.  It  is  obvious  that  there  was  simply  no  fit  of  any  of 
the  state  vector.  It  should  be  noted  that  the  plot  of 
average  horizontal  stabilizer  position  although  appearing  to 
be  very  rough  was  a  result  of  the  plot  scale  only  and  was  in 
essense  relatively  smooth. 

The  parameter  identification  results  using  the  5  SPS 
data  indicated  it  was  possible  to  identify  the  phugoid  mode 
parameters.  Referring  to  Table  VIII,  it  can  be  seen  that  the 
estimated  derivatives  using  the  5  SPS  data  converged  to 
appropriate  (sign  and  magnitude)  stability  derivative 
values. 


Derivatives    Z       X       X       X       M       X A 

u       oc  u       g       u       ^e 

Estimate    -.0002       0     -.006    -45.5   -9.6x10-s   0 

Standard    1.2x10~6   fixed    .0002    fixed   5.x10~6    fixed 
Deviation 

Parameter   -8.7x10-6  fixed   -.0016    fixed  -.00002     fixed 
Step  Size 

Eigenvalue  1.82  x  10*2 

Spread 

Likelihood 

Function  -  -1.49  x  10* 

Value 

Time 

History        0~  =    0.00267      o~     =    5.53    <r    =  0.00343 

Matches         oe  u  g 

<T     =  0.014  CT    =    0.92 
e  n 

z 

Table  VIII.   SCIDNT  Estimates  of  Longitudinal  Phugoid 
Stability  Derivatives  using  5  SPS  Data 


The  standard  deviations  of  the  parameter  estimates  were  much 
smaller  than  the  parameter,  and   the   parameter   step   sizes 
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indicated  that  the  estimates  had  converged  to  the  "best" 
estimates  obtainable  with  the  data.  It  should  be  noted  that 
standard  deviations  for  u  and  6  were  5.53  ft  sec-1  and  0.014 
rad  respectively,  a  significant  reduction  in  standard 
deviation  from  those  obtained  with  2  SPS  data 

The  time  history  plots  using  5  SPS  data  showed 
relatively  close  agreement  between  the  airplane  and  model 
time  histories.  There  was,  however,  some  model  instability 
noted  during  the  first  30  seconds  of  the  time  histories.  It 
was  concluded  that  these  instabilities  were  due  to  the 
relatively  slow  sampling  rate  during  the  initial  control 
input. 

Prom   SCIDNT   the   phugoid   damping   ratio  (  f  )  and 

freguency    (ct>   (rad  sec-1))    were    0.0486    and    0.06 

n 

respectfully.    Available  F-14  data  was  obtained  in  the  auto 

sweep  mode  with  SAS  on  while  SCIDNT  data  was  obtained  in 
auto  sweep,  SAS  off.  It  was  therefore  not  possible  to  make 
any  meaningful  comparisons.  However,  it  was  noted  that 
between  10,000  ft  and  35,000  ft  the  phugoid  period  varied 
between  106  sec  and  72  sec.  The  SCIDNT  freguency  resulted  in 
a  period  of  92.53  sec  which  when  considered  in  light  of 
available  NPE  data  appeared  reasonable. 

It  was  concluded  that,  within  the  scope  of  this 
research,  the  paramater  identification  program  SCIDNT 
provided  a  reasonable  estimate  of  the  phugoid  longitudinal 
parameters. 


C.   SIMULTANEOUS  SHORT  PERIOD  AND  PHUGOID  PARAMETER 
IDENTIFICATION 

1 .   Airplane  T es t  Conditions 

Initial  trimmed  test  conditions  for  the  short  period 
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and  phugoid  maneuvers  are  presented  in  Table  IX.  Prior  to 
the  control  input  the  airplane  was  trimmed  for  steady  state 
"hands  off"  flight  [5].  Since  ths  short  period  input  was 
applied  after  the  phugoid  input,  the  two  test  conditions  are 
different. 


Altitude     MACH       CG       Loading 
(Ft  HPC)  (%  MAC) 

Phugoid  Input        30  326      .6217     10.86         B 

Short  Period  Input   30  703      .5963     10.86         B 

Table  IX.   Short  Period  and  Phugoid  Test  Conditions 


2 .   Method  Of  Investigation 

The   airplane   model  equations  used  were  those  shown 

in  Appendix  B.  The  derivatives  X   ,  X  -  ,    and  X    were   fixed 

oc    °e        q 

as  described  in  the  previous  section. 

The  maneuver  investigated  in  this  section  was 
designed  specifically  to  provide  the  means  of  identifying 
the  short  period  and  phugoid  parameters  from  one  maneuver 
or,  if  you  will,  simultaneously.  Initiation  of  the  phugoid 
maneuver  was  accomplished  in  accordance  with  standard  test 
techniques,  [5],  Upon  establishing  the  airplane  in  the 
phugoid  maneuver,  a  doublet  input  was  applied  to  excite  the 
short  period  mode.  The  airplane  was  then  allowed  to  fly 
through  at  least  one  phugoid  cycle. 

Identification  of  the  short  period  parameters  was 
essentially  routine  and  straightforward.  The  proceedures 
used  in  estimating  the  parameters  were  those  described  in 
the  previous  short  period  section.  The  computer  runs  were 
made  with  the  number  of  iterations  and  step  cuts  both 
initialized  at  ten.  Although  this  short  period  input  was 
superimposed   on  the  phugoid,  convergence  was  obtained  early 
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in  program  execution.  Iterations  conducted  under  these 
conditions  converged  very  quickly  and  executed  step  cuts 
until  the  allowable  number  was  exceeded.  Likelihood  function 
change  during  the  short  period  iterations  was  similar  to 
that  shown  in  Figure  2. 

Based  on  the  results  of  the  previous  section  it  was 
concluded  that  data  recorded  at  5  SPS  would  most  likely 
yield  reasonable  phugoid  parameter  estimates.  The  short 
period  stability  derivatives  obtained  from  this  run  were 
fixed  as  constants  during  the  phugoid  iterations.  Computer 
runs  were  then  made  with  the  number  of  iterations  and  step 
cuts  both  initialized  at  ten. 

Figure  6  presents  the  likelihood  function  value 
during  iteration  of  the  phugoid  stability  derivatives. 
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Figure  6.   Likelihood  Function  During  Iteration 

For  Phugoid 
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The  likelihood  function   change   for   the   phugoid   maneuver 

displayed     some     characteristics     heretofore    not 

encountered.  The  characteristics  shown  are   typical   of   the 

results   obtained   with   the   available   data.  At   some  time 

during  the  iteration  sequence  the  likelihood  function   would 

exhibit   a   jump   discontinuity,  in   this   case   it  occurred 

between  the  first  and  second   iterations.  Program   iteration 

then   appeared  to  proceed  in  a  normal  manner  until  the  tenth 

iteration  at  which  time  the  program  terminated.  From   Figure 

6   it   was  obvious  that  convergence  had  not  been  obtained  at 

program  termination. 

In  an  attempt  to  correct  the  iteration  instabilities 

the   phugoid   stability   derivatives   Z  ,  X  ,  and   M     were 

u    u        u 

estimated   simultaneously   with  Z  ,    X  ,  H  ,  and  6  .  Figure  7 

o   o   o       o 

presents  the  likelihood  function  during  a  typical   iteration 
for  seven  parameters. 
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Figure  7.   Likelihood  Function  During  Iteration  For 

Phugoid  (Stability  Derivatives,  Z  ,  X  ,    M  ,  and  6  ) 

o    o   o       o 


It  can  be  seen  that  the  likelihood  function  value  converged 
rapidly  to  a  steady  state  value  with  no  apparent 
discontinuities  or  nonlinearaties. 

3 •   Data  Analysis 

Data  from  a  typical  short  period  run  superimposed  on 
the  phugoid  are  presented  in  Table  X.  A  typical  time  history 
fit  is  presented  in  Appendix  D,  Figure  6. 
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Derivatives    Z       Z       M       H       Z^      Mx, 

cc       g       oc      g       ffe      "e 

Estimate      -.4285   1.048    -1.54    -.665    0.19     4.464 

Standard      0.0049   0.0074   0.014    0.011    0.0082   0.0312 
Deviation 

Parameter     0.0149   0.0119   -.0166   0.0355   0.0173   -.1815 
Step  Size 

Eigenvalue  1.25  x  105 

Spread 

Likelihood 

Function  -4.36  x  103 

Value 

Time 

History        (T  =  0.00188      Q~     =  0.8822  <T    =  0.0036 

Hatches         cC  u  g 

(T  =  0.0048  (T     =    1.82 
6  n 

z 

Table  X.   SCIDNT  Estimates  of  Longitudinal  Short  Period 

Stability  Derivatives 


The  short  period  parameter  estimates  displayed 
characteristics  similar  to  those  observed  when  estimating 
the  short  period  parameters  from  a  short  period  input 
only.  The  time  history  plots  showed  good  agreement  between 
the  model  and  actual  airplane  data.  It  was  noted  that 
airspeed  was  not  constant  since  the  short  period  was 
superimposed  on  the  phugoid.  This  varying  airspeed  did  not 
however,  have  any  noticable  effect  on  the  determination  of 
the  short  period  parameters. 

The  results  for  the  phugoid  portion  of  the  maneuver 
were  at  times  in  disagreement  with  theory.  From  theory  it 
can  be  shown  that  the  speed  damping  derivative,  which  acts 
along  the  negative  X  axis  increases  with  speed  and  is 
therefore  negative  [ 6 ]•  During  the  iteration  to  determine 
the   phugoid   stability   derivatives   it   was  noted  that  the 
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speed   damping   derivative   (X  ) ,  became   positive   on    the 

u 

seventh    iteration   and  remained   positive   until   program 
termination,  Table  XI. 


Derivatives    Z       X       X       X       H       X~ 

u       OC       u       g        u        °e 

Estimate    -.0002     0        .005    -43.85  -7.8x10-*   0 

Standard    1.1x10-6   fixed    .0002    fixed   3.5x10~6   fixed 
Deviation 

Parameter   5.38x10~6  fixed   4.8x10-5  fixed  1.2x10~s    fixed 
Step  Size 

Eigenvalue  6.2  x  1012 

Spread 

Likelihood 

Function  -1.005  x  10* 

Value 

Time 

History        CT  =  0.004        <T      =    1.97   GT     =  0.0047 

Matches         oc  u  g 

(T   =  0.031   CT  =  1.48 
6  n 

z 

Table  XI.   SCIDNT  Estimates  of  Phugoid  Stability 

Derivatives  (Positive  X  ) 

u 


The    time    history    plots    were     in     reasonable 

agreement,  Appendix  D  Figure  7.  The  positive  value  of  X   was 

u 

not  however  considered  reasonable   since   this   value   would 

imply   decreasing  drag  with  increasing  airspeed.  The  X   term 

u 

was  however  negative  through  the  sixth  program  iteration.  As 

mentioned    previously   the    likelihood    function   during 
iteration   showed   definite   signs   of   instability,  it   was 
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therefore  suspected  that  the  speed  damping  sign  change  was 
due  to  this  program  instability.  The  program  was 
therefore,  run  again  and  terminated  at  the  sixth  iteration 
thereby  producing  a  negative  speed  damping  term,  Table  XII. 


Derivatives    Z       X       X       X       M       Xrf 

u       cc       u       g       u       ve 

Estimate    -.00008      0      -.00032  -43.85  0.00017      0 

Standard    2.6x10~6   fixed    0.00047  fixed   6.7x10^*   fixed 
Deviation 

Parameter   4.0x10~s   fixed   -.0059    fixed  0.000058    fixed 
Step  Size 

Eigenvalue  8.53  x  10*1 

Spread 

Likelihood 

Function  -7.68  x  103 

Value 

Time 

History        0~  =    0.0071       Q-     -    3.94    <T    =0.0053 

Matches         oc  u  g 

C     =    0.0366  CT    =    2. 14 
6  n 

z 

Table  XII.   SCIDNT  Estimates  of  Phugoid  Stability 

Derivatives  (Negative  X  ) 

u 


With  the  iterations  thus  restricted  to  six  the  value  for   X 

u 

was   maintained  at  a  negative  value.  Time  history  plots  for 

the  negative  X  case  are   presented   in   Appendix   D,  Figure 
u 

8.  It  was  noted  that  there  was  no  appreciable  improvement  in 

the  time  history  jlots.  In  fact  from  Tables  XI   and   XII   it 
can   be   seen  that  the  standard  deviation  of  the  state  error 
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vectors   was   degraded   for   the   negative    X    case.  This 

u 

condition   in  fact  is  however,  not  very  startling  when  it  is 

realized  that  the  program  had  demonstrated  instabilities  and 

had  not  converged. 

In  an  effort  to  reconcile  these  inconsistencies  and 

lack   of   convergence   the   three   phugoid   parameters   were 

estimated   along   with   Z  ,  X  ,    h    ,    and   6  .  This   technique 

o    o   o        o 

yielded  reasonable  stability  derivatives  and  acceptable  time 

history  plots.  Stability  derivative  data  provided  by  the 
seven  parameter  identification  are  shown  in  Table  XIII,  time 
history  plots  are  presented  in  Appendix  D,  Figure  9. 


Derivatives     Z        X       X        X        M        X^ 

u       oc       u       g       u       °e 

Estimate    -.00019      0      -.0029   -43.85  -.000141     0 

Standard    1.4x10~6   fixed    .00026   fixed   5.2x10~&   fixed 
Deviation 

Parameter   1.23x10~7  fixed   6.3x10~s  fixed  -4.02x10~6  fixed 
Step  Size 

Eigenvalue  3.7  x  1013 

Spread 

Likelihood 

Function  -4.23  x  10^ 

Value 

Time 

History        CT  =  0.0045       0~     =  4.17    <r    =0.00523 

Matches         cc  u  g 

<r    =0.031  <r    =  1.86 
e  n 

z 

Table  XIII.   SCIDNT  Estimates  of  Phugoid  Stability 

Derivatives  (Parameters  Estimated  With  Z  ,  X  ,  M  ,  and  e  ) 

o   o   o       o 
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The  parameter  standard  deviations  were  small  and 
within  engineering  tolerance.  The  parameter  step  sizes  were 
quite  small  and  showed  further  evidence  of  solution 
convergence.  The  standard  deviation  of  the  calculated  model 
response  with  respect  to  the  measured  airplane  response  was 
acceptable.  The  time  history  plots  provided  further 
indication  of  parameter  estimate  validity. 

The  problems  associated  with  the  phugoid  using  the 
simultaneous  technique  could  have  been  a  result  of  program 
deficiencies  or  airplane  anomalies  or  a  comDination  of 
both.  The  flight  conditions  at  which  the  data  was  obtained 
were  areas  in  which  the  phugoid  mode  was  lightly  damped  or 
neutrally  damped.  In  addition  the  airplane  wing  was  set  to 
the  auto  sweep  mode.  It  was  felt  that  either  of  these 
conditions  could  have  caused  the  program  instabilities 
noted.  Attempts  to  obtain  new  data  were  not  possible  due  to 
airplane  non-availability  as  mentioned  earlier. 

It  was  concluded  that  extraction  of  the  short  period 

stability  parameters  in  the   presence   of   the   phugoid   was 

easy,  straight  forward,  and  yielded  results  similar  to  those 

obtained  from  pure   short   period   data.  Estimation   of   the 

phugoid   stability   parameters   was   possible  when  they  were 

estimated  in  conjunction  with  Z  ,    X  ,  H  ,  and  9  . 

o    o   o       o 

It   was   recommended   that   a   new   set   of   data  be 

obtained  at  a  flight  condition  where  the  phugoid  is  at  least 

moderately   damped  and  that  this  data  be  analyzed  to  resolve 

the  present  anomalies. 
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Ill-   CONCLUSIONS 


A.   GENERAL 

The  maximum  likelihood  identification  program  SCIDNT 
showed  great  promise  in  obtaining  the  short  period  and 
phugoid  stability  derivatives  from  one  maneuver, 
"simultaneously". 


B.   SPECIFIC 

Within  the  scope  of  this  research,  the  parameter 
identification  program  SCIDNT  was  easy  and  straight-f crward 
to  use  in  determining  the  short  period  mode. 

Meaningful  parameter  estimates  were  obtainable  only  when 
started  from  an  initial  steady  state  condition. 

Within  the  scope  of  this  research,  the  parameter 
identification  program  SCIDNT  provided  a  reasonable  estimate 
of  the  phugoid  longitudinal  parameters. 

Extraction  of  the  short  period  stability  parameters  in 
the  presence  of  the  phugoid  was  easy,  straight-forward,  and 
yielded  results  similar  to  those  obtained  from  pure  short 
period  data. 

Estimation  of  the  phugoid  stability  parameters  was 
possible   when   they   were   estimated   in   conjunction   with 

Z  ,  X  ,  H  ,  and  9  . 
o   o   o       o 


72 


IV.   RECOMMENDATIONS 

It  is  recommended  that  a  new  set  of  data  be  obtained  at 
a  flight  condition  where  the  phugoid  is  at  least  moderately 
damped  and  that  these  data  be  analyzed  to  resolve  the 
present  anonalies. 
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APPENDIX  A 


TEST  AIRCRAFT  DESCRIPTION 


The  F-14  airplane  is  a  supersonic,  two-place, 
twin-engine,  swing-wing,  air-superiority  fighter  designed 
and  manufactured  by  Grumman  Aerospace  Corporation.  In 
addition  to  its  primary  fighter  role,  carrying  missiles  and 
an  internal  20-miilimeter  gun,  the  airplane  is  designed  for 
fleet  air  defense  and  ground  attack  missions.  Photographs  of 
the  test  airplane  are  shown  in  Figures  1  and  2.  A  three-view 
drawing  is  presented  in  Figure  3. 

The  forward  fuselage,  containing  the  crew  and  electronic 
equipment,  projects  forward  from  the  mid  fuselage  and  wing 
glove.  Outboard  pivots  in  the  highly  swept  wing  glove 
support  the  movable  wing  panels,  which  incorporate  integral 
fuel  cells  and  full-span  leading  edge  slats  and  trailing 
edge  flaps  for  supplemental  lift  control.  In  flight,  the 
wings  may  be  varied  in  sweep,  area,  camber,  and  aspect  ratio 
by  selection  of  any  leading-edge  sweep  angle  between  20°  and 
68°.  Wing  sweep  can  be  automatically  or  manually  controlled 
to  optimize  performance  and  thereby  enhance  aircraft 
versatility.  Separate  variable-geometry  air  inlets,  offset 
from  the  fuselage  in  the  glove,  direct  primary  airflow  to 
two  TF30-P-412A,  dual-axial-compressor,  turbo-fan  engines 
equipped  with  afterburners  for  thrust  augmentation.  The 
displaced  engine  nacelles  extend  rearward  to  the  tail 
section,  supporting  the  twin  vertical  tails,  horizontal 
tails,  and  ventral  fins.  The  middle  and  aft  fuselage,  which 
contains  the  main  fuel  cells,  tapers  off  in  depth  to  the 
rear  where  it  accommodates  the  speed  brake  surface  and 
arresting  hook.  Retractable  vanes  in  the  glove  leading  edges 
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extend  to  supplement  lift  and  compensate  for  changes  in  the 
aircraft  aerodynamic  center.  All  control  surfaces  are 
positioned  by  irreversible  hydraulic  actuators  to  provide 
desired  control  effectiveness  throughout  the  flight 
envelope.  Stability  augmentation  features  in  the  flight 
control  system  enhance  flight  characteristics  and  provide  a 
more  stable  and  maneuverable  weapons  delivery  platform.  The 
tricycle-type,  forward-retracting  landing  gear  is  designed 
for  nose  gear  catapult  launch  and  carrier  landings.  Missiles 
and  external  stores  are  carried  from  eight  hardpoint 
stations  on  the  center  fuselage  between  the  nacelles  and 
under  the  nacelles  and  wing  glove;  no  stores  are  carried  on 
the  movable  portion  of  the  wing.  The  fuel  system 
incorporates  both  inflight  and  single-point  ground  refueling 
capabilities. 

A   more   complete   description   of  the  F-14A  airplane  is 
given  in  Ref.  8. 
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3/4  FRONT  VIEW 


FRONT  VIEW 


F-14A  AIRPLANE 


Figure  1 
APPENDIX  A 
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SIDE  VIEW 


F-14A  AIRPLANE 


Figure  2 
APPENDIX  A 
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THREE    VIEW    DRAWING 
F-14A    AIRPLANE 


Figure   3 
APPENDIX    A 
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APPENDIX  B 
DERIVATION  OF  THE  AIRPLANE  MODEL  AND  MEASUREMENT  EQUATIONS 

A  detailed  derivation  of  the  complete  aircraft  equations 
of  motion  is  shown  in  Ref.  6.  A  summary  of  the  assumptions 
made  during  that  derivation  is  shown  below  for 
completeness.  These  assumptions  lead  to  a  set  of  force  and 
moment  equations. 

*   Summary  of  Derivation  Assumptions   * 

1.  The  airframe  is  a  rigid  body 

2.  The  earth  is  fixed  in  space 

3.  The  earth's  atmosphere  is  fixed  with  respect  to  the 
earth 

4.  The  mass  of  the  airplane  remained  constant  for  the 
duration  of  any  particular  dynamic  analysis 

5.  The  xz  plane  is  a  plane  of  symmetry 

6.  The  disturbances  from  the  steady  flight  condition  is 
small  enough  so  that  the  products  and  squares  of  the  changes 
in  velocities  are  negligible  in  comparison  with  the  changes 
themselves 

7.  The  disturbance  angles  are  small  enough  such  that  the 
sines  of  these  angles  may  be  set  equal  to  the  angles  and  the 
cosines  set  equal  to  one.  Products  of  these  angles  are  also 
approximately  zero  and  can  be  neglected 

8.  Since  the  disturbances  were  small,  the  change  in  air 
density  encountered  by  the  airplane  during  any  disturbance 
can  be  considered  zero 
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9.  During  the   steady   flight   condition,  the  airplane   is 

flying   with  wings  level  and  with  all  components  of  velocity 

zero  except  U   and  W 
o      o 

10.  Body  axis  system  is  assumed 


Dividing  the  force  equations  by  the  mass  m  and  the 
moment  equations  by  the  appropriate  moments  of  inertia 
yields  terms  of  the  form: 


J_lX_u_        and   1   3L  r.        (  B1  ) 
m  7u  T*  fir 

XX 

Replacing   the  first   term   by  X   and  the  second  term  by  L 

u  r 

simplifies   the   notation.  These    quantities   are    called 

"dimensional   stability   derivatives".  Using  this  short  hand 

notation   the   equations   of   motion    applicable    to   the 

longitudinal   problem  can  be  written  as  follows:  (a  complete 

set  of  eguations  is  presented  in  Ref.5) 


U  +  W  q  +  ge  COS  0   =Xu+X,U+Xq+Xq 
o  o     u      u      q      q 

♦  Xw+X.w  +  X^cT   +  X^  O       +  Xy  O 
w     w      ee     ee     ee 


F   F     °F   F     °B   B      B   B 


Xu      (j      +    cos    £   T    u    +    cos  £    Tx  tf 

°B      B  ^      u  o       °RPM      RPM 


(    B2    ) 
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w  -   U   g   +   g6   sin  e      =   Z   u  +    Zru  +   Z   g   +   Z.g 

o                              o           u             u  g             q 

♦Zw+Z.v+Z.     0      ♦   Z'6     <5      +  Z^    *<5 

w             w               e      e             ee  ee 


+   Z.5     (5      +    Zi     <5      +   Ztf     6      +    Z£     <S 


F      F  ~F      F  B      B  B       B 

+    Z5f     #      -    (sin  £  )T    u   -    (sin  £    )  T^         ^5  (    B3    ) 

UR         R  O  n  O  ^PPMRPM 


B       B  °  U  °  RPM       RPM 


g=Mu+M.u+Mg+M,g+Mw+Mtw+M£     tf 
u  u  g  g  w  w  ee 

+    Mi      6      +    MS      0      +    Mtf     <J      +    M>     £      ♦    M^     (5 
°e      e  °e      e  °F      F  °F      F  °B      B 

z   m  z   m 

..  J  J 

+    M^       6       +    MV     a     + T    u    + T  ,*  (S  (    BU    ) 

°BB  °BBI  u  I  °  RPM       RPM 

yy  yy 


In  order  to  simplify  these  eguations  it  was  assumed  that 
all  disturbances  were  small  such  that  the  higher  order  terras 
were  negligible.  Also  it  was  assumed  that  the  flaps  and 
speedbrakes  were  retracted  and  that  thrust  effects  were 
negligible.  With  these  further  assumptions  the  w  eguation 
became: 


w   -    0    g    +    ge    sine      =    Zu+Zg    +    Zw    +    Zx#  (B5) 

o  o  u  g  w  e      e 


It  was  further  assumed  that  the  flow  was 
guasisteady.  Because  of  this  all  derivatives  with  respect  to 
the  rates  of  changes  of  velocities  were  omitted  with  the 
exception  of  those  with  respect  to  w  which  were  retained  to 
account  for  the  effect  on  the  horizontal  tail  of  the 
downwash  from  the  wing.  If  a  wing  moving  through  an  air  mass 
with   an  angle  of  attack cC   and  a  steady  forward  velocity  U 
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is  suddenly  given  a  vertical  velocity  change  w  without 
changing  the  forward  velocity,  the  angle  of  attack 
perturbation  is  as  shown  in  Figure  1. 


(Perturbation  in  w) 


Figure  1.   Angle  of  Attack  Due  to  Vertical  Velocity 

Perturbation 


Thus   changes  in   angle  of  attack,     ,  are  proportional  to 
the  vertical  velocity  perturbation  w.  Thus  we  can  write, 


w  +  W 
o 
tan(0C   +  °c  )  =  ,  , 

o        Yu*2  +  (w  +  W  )  2 


(  B6  ) 


and  by  the  small  angle  assumption 

W 
o 
and  oc     =  


OC  = 


D 


D 


(  B7  ) 


Taking  the  derivative  yields, 


w  =  U  OC 
o 


(  B8  ) 


also  it  can  be  shown  that, 


3Z  aZ 

=  Limit  

j}w    w — ►O   w 


(  B9  ) 
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and    substituting   eguation    B    8    yields, 


3Z  1      aZ 

-5—=   Limit 

o    U       « 
o 


aw     oc- 


u 


CO 


(    B10    ) 


also   if    we   redefine   the   Z    dimensional    derivatives   such   that, 
typically, 


1  ,1  .1 

Z        =  —  z  ,      z     =  Z  ,      Z     =  —  z 

ccu<x  uuu  g      u     q 

o  o 


U         *e 


(    B11    ) 


we   can   now   substitute   into  equation  B  3  which  yields  the 
equation, 


CC  =  z'  OC   ♦  Z'  u  +  Z'  q 
cC        u      q 


9  / 

—  e  sm  e  +  z*     a 

0         o    °e   e 


(  B12  ) 


where  Z   =  1  +  Z 
<3       g 


A  similar   expression   was   obtained   for   u   using  the 
simplifying  assumptions  and  the  expressions, 


3X          1   aX     1 
=  Limit =  X 


dw    oc— *o  U   ^    U    oc 
o        o 


(  313  ) 


and  the  expression  for  u  was  therefore  written, 


U=X      CC+Xu+Xq   -g0    cos6      +    X*      6 
<X  u  q  o  e      e 


(    B14    ) 
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where  X'  =  (X   -  cc  U  ) 
q      g      o  o 

Finally  the  expression  for  q  followed  directly  as, 
q  =  M  OC     +M  U  +  ;1  q  +M«  6 


oc 


(  B15  ) 


e   e 


The  model  equations  for  the  longitudinal  case,  including 
an  angle  of  attack  due  to  gusts  term,  were  written  in  matrix 
form  as, 


x  =  Fx  +  Gu  +  Pw 


(  B16  ) 


where 


state 
vector 


x  = 


OC 
u 

g 

o 

oc 

g 


(angle  of  attack  perturbation) 
(longitudinal  velocity  perturbation) 
(pitch  rate  perturbation) 
(pitch  attitude  perturbation) 
(angle  of  attack  due  to  gusts) 


control 
vector 


u  =  I  6"     (horizontal  elevator  deflection  froi 
L  e  J    trim) 


F  = 


Z* 


oc 


0 
0 


z' 

q 

x' 

q 

q 
1 

o 


g/U   sine 
o     o 

-g  cos  6 


0 
0 


cC 


oc 


oc 


-CO 

c 


u 


Note:   z'   in   the   model   is   equivalent  to  the  equation  of 

q 

A 

motion  term  (1  +  Z  ) 

q 


8U 


g  = 


L   J 


W  = 


The  measurement  equations  written  in  matrix  form  were, 


y  =  Hx  +  Du  ♦  v 


(  317  ) 


HX  +  Du  = 


K   (  cC  +  CC   )  + 

*      g 

u 


g 
e 


V  *c* 


-0  z   (oc  +  ocjz^  6 
o   °C         g    ve   e 


v  = 


n 


n 


Note:  Selected  stability  derivatives  are  primed  since 
they  differ  from  the  notation  used  in  Ref.  6,  however  they 
are  exactly  equivalent  to  the  derivatives  used  in  SCIDNT. 
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APPENDIX  C 


AIRPLANE  AND  GROUND  DATA  SYSTEMS 


A.   AIRPLANE  DATA  SYSTEM 


1.   General 


The  demonstration  F-14A  ship  number  8  BU.  NO.  157937 
was  the  test  airplane  utilized.  A  detailed  description  of 
the  airplane  instrumentation  system  is  contained  in  Refs.  7 
and  9.  The  airplane  was  modified  to  be  compatible  with  the 
Naval  Air  Test  Center  Real-Time  Telemetry  Processing  System 
(RTPS) .  The  instrumentation  systems  and  measurements 
selected  for  use  in  the  test  airplane  were  primarily 
designed  to  satisfy  the  Flight  Test  requirements  fcr  the 
F-14A  Performance  Evaluation  and  Demonstration.  The  test 
airplane  was  not  entirely  a  production  configured.  F-14A  in 
that  some  miscellaneous  cockpit  components  not  required  for 
flight  testing  were  removed.  The  production  nose  radome  was 
replaced  by  a  fiberglass  and  metal  radome  fitted  for  an 
instrumentation  noseboom  assembly.  Additionally,  some 
production  equipment  was  reworked  in  order  to  accomodate 
various  instrumentation  subassemblies. 

2-   Nose  Boom  Instrumentation 

The   nose  boom  assembly  included  Angle  of  Attack  and 
Angle  of  Side  Slip  Vanes  and  a  Pitot  /  Static  Head. 
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3«   Cockpit  Instrument at ion 

The  Forward  Cockpit  instrumentation  was  primarily 
pilot  oriented  and  included;  Angle  of  Attack,  and  Sideslip 
Indicators,  Airspeed,  Altimeter,  Machmeter,  Center  of 
Gravity  Normal  Acceleration  and  an  Instrumented  Control 
Stick.  The  Aft  Cockpit  Instrumentation  included  left  and 
right  engine  RPM  (N1),  Main  Fuel  Flow,  A/B  Fuel 
Flow,  Outside  Air  Temperature  and  the  Instrumentation  Master 
Control  Console. 

4 .  Compartment  Instrumentation 

The  instrumentation  subcomponents  were  mounted 
throughout  the  airplane  in  various  available  locations 
including,  the  Nose,  Forward  Compartments,  and  Equipment 
Compartments.  The  bulk  of  the  instrumentation  equipment  was 
located  in  the  forward  equipment  compartments. 

5.  System  Description 

The  instrumentation  system  consisted  of  the 
following  subsystems: 

o  PCM  Signal  Conditioning 

o  PCM 

o  Timing 

o  Magnetic  Tape  Recording 

o  Telemetry 

o  Power  and  Distribution 

o  Ground  Support  Equipment 

The  system  had  the  capacity  to  record  on  magnetic 
tape  and  to  telemeter  to  the  ground-data-reduction-ccmplex 
all  the  PCM  data  for  real-time  analysis.  The  PCM  Signal 
Conditioning  Subsystem  consisted  of  active  pre-sample  filter 
chassis  and  seme  special  conditioning  required  for  certain 
measurements.  The  pre-sample  filters  had  cut-off  frequencies 
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of  5,  10,  20,  40,  and  80  Hz,  The  PCM  subsystem  consisted  of 
an  analog  multiplexer,  analog  to  digital  converter,  digital 
multiplexer,  and  the  required  control  logic.  The  subsystem 
had  the  capacity  of  accepting  up  to  332  analog  parameters  at 
frequency  response  from  dc  to  40  Hz.  Discrete  data  (On/Off 
signals)  and  a  limited  amount  of  digital  data  and  time  words 
were  time-division  multiplexed  into  the  output  data 
stream.  The  output  data  was  presented  as  a  succession  of 
twelve-bit  (including  parity)  parallel  data  words  in  NRZM 
format  which  were  recorded  on  the  digital  tape.  A  serial 
output  in  NRZL  format  for  entry  into  the  "L"  Band  telemetry 
subsystem  was  provided.  The  PCM  format  is  contained  in  Ref.- 
10.  An  IRIG-E  Time  Code  Generator  was  installed  to  provide  a 
common  time  correlation  for  magnetic  tape  and  telemetry 
data.  The  Magnetic  Tape  Recording  Subsystem  consisted  of  a 
PCM  Recorder  which  had  the  capability  to  record  the  PCM  data 
and  Time  Code  information.  The  Digital  Recorder  contained  a 
sixteen  tracks  IRIG  compatible  head  capable  of  storing  one 
hour  of  digital  data  on  one  inch  tape  at  a  tape  speed  of  30 
IPS.  The  data  was  recorded  in  parallel  using  saturation 
recording  techniques.  This  portion  was  capable  of  bit 
densities  of  1000  bits  per  inch.  The  Telemetry  Subsystem 
consisted  of  a  signal  combining  unit,  an  "L"  Band 
Transmitter  and  two  blade  antenna.  The  operating  frequencies 
could  be  tuned  to  coinside  with  those  in  use  at  various  test 
facilities.  The  entire  Flight  Test  Instrumentation  System 
was  powered  from  the  aircraft  primary  (generating) 
system.  Should  either  the  left  or  right  engine  or  generator 
fail,  the  instrumentation  system  would  continue  until  such 
time  as  both  generators  failed.  The  distribution  system 
consisted  of  a  28V  dc  transformer/rectifier  (T/R)  ,  and  a 
power/control  distribution  box,  which  provided  electrical 
system  protection  and  isolation  from  the  aircraft  power 
systems  through  individual  circuit  breakers.  The  equipment 
required  to  perform  checks  on  the  installed  instrumentation 
systems   was   configured   in   mobile  test  carts.  These  carts 
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required  only  60  cycle  power  for  operation  and  provided  the 
operator  with  a  tool  for  conducting  pre-flight  and  post 
flight  tests  on  the  tape-  recording  system. 


B.   GROUND  DATA  SYSTEM 

The  Naval  Air  Test  Center,  Patuxent  River,  Maryland, 
employs  a  xieal-Time  Telemetry  Processing  System  (RTPS)  as 
universal  support  for  testing  on  several  types  of  aircraft 
and  weapons  systems, Ref.  10.  The  computer  system,  designed 
by  Xerox,  consists  of  a  Sigma  9  general  purpose  computer  as 
the  central  processor  supported  by  three  Sigma  3  computers 
performing  the  data  input  and  output  processing.  RTPS 
provides  the  capability  to  process  two  complete  and 
independent  tests  (data  "streams")  simultaneously  from 
either  tape  or  telemetry  sources.  These  "streams"  can 
support  totally  unrelated  and  asynchronous  test  flights  or 
can  support  setup  for  one  flight  while  processing  another 
flight. 

Each  Data  Channel  accepts  up  to  512  channels  of  PCM  data 
and  32  parameters  of  FM  data.  This  data  is  processed 
(individually  or  mixed)  by  a  Sigma  3  and  preprocessing 
Programmed  Algorithm  Unit  (PAU) .  Up  to  50,000  values  of 
Engineering  Units  (EU)  data  per  second  can  be  fed  to  the 
Sigma  9  while  simultaneously  being  stored  on  digital  tapes 
for  later  playback.  Two  large  screen  interactive  graphic 
display  stations,  that  provide  real-time  "hands  on"  data 
analysis  and  presentation,  are  controlled  by  the  Sigma  3 
Display  Processor.  Batch  jobs  can  be  processed  when  the 
real-time  and  non  real-time  operational  modes  are  not 
utilizing  the  resources  of  the  system  to  100  percent.  The 
batch  job  is  defered  if  additional  demands  are  made  by 
either  non  real-time  or  real-time  tasks  during  batch  job 
processing. 

The   specific  flight  test  data  utilized  by  the  parameter 
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identif icaticn   program   were,  V  ,  ALPHA   FRL,  pitch    rate, 

T 

pitch   attitude,   NZ   NET,  and  average  horizontal  stabilizer 

position.  The  raw  data  (digital)  was  used  to  generate  an  EU 
tape  for  input  to  the  parameter  identification  program,  in 
generating  the  EU  tape  corrections  were  made  to  the  measured 
values  as  required  in  accordance  with  Ref  11.  For  example 
corrections  tc  angle  of  attack  included  nose  boom  bending 
moments,  fuselage  bending  moments  and  pitch. 
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